Prueba Claro

Frank Montenegro - frandak2@gmail.com

importamos las librerias que vamos a usar

import pandas as pd
import geopandas as gpd
import libpysal as lp
import matplotlib.pyplot as plt
import seaborn as sns
import geoplot as gplt
import geoplot.crs as gcrs
import plotly.express as px
import rasterio as rio
import numpy as np
import contextily as ctx
import shapely.geometry as geom
import folium
import folium.plugins as plugins
import os
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')# Configuración warnings

redireccionamos la ruta del notebook

# import os
g = os.getcwd()
g
path = r'D:\repos\GEO_CLARO'
os.chdir(path)

Procedemos a leer los archivos shapefile

de antemanos sabemos que el shape MGN es la capa de centros poblados del igac, en este shape tenemos imforacion sobre la forma, tamanio y area de cada centro poblado del pais

points_late = gpd.read_file("data\points_late.shp")
points_late.head(5)

	FECHA	TEST_ID	N_OPERADOR	DOWNLOAD_M	UPLOAD_MBP	CLASS_LATE	localidad	MES	DIA	geometry
0	2021-01-28	6978771403	4	10.410156	0.640625	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.09600 4.58700)
1	2021-01-19	6945885369	4	10.294922	4.988281	0	ANTONIO NARIÑO	January	Tuesday	POINT (-74.12100 4.59200)
2	2021-01-08	6907797685	4	28.707031	4.711914	0	ANTONIO NARIÑO	January	Friday	POINT (-74.10000 4.59000)
3	2021-01-21	6953823661	4	53.502930	7.201172	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.09500 4.58500)
4	2021-01-21	6952553476	4	26.608398	6.316406	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.12100 4.59300)

Que podemos averiguar con este dataset

- crear un mapa del total de comparendos mensuales por localidad?

- crear un densidad de comparendos anual?

points_late['FECHA'] = pd.to_datetime(points_late['FECHA'])## transformamos la columna fecha
points_late['MES'] = pd.Categorical(points_late['MES'], categories=['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September'])
points_late.head(5)

	FECHA	TEST_ID	N_OPERADOR	DOWNLOAD_M	UPLOAD_MBP	CLASS_LATE	localidad	MES	DIA	geometry
0	2021-01-28	6978771403	4	10.410156	0.640625	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.09600 4.58700)
1	2021-01-19	6945885369	4	10.294922	4.988281	0	ANTONIO NARIÑO	January	Tuesday	POINT (-74.12100 4.59200)
2	2021-01-08	6907797685	4	28.707031	4.711914	0	ANTONIO NARIÑO	January	Friday	POINT (-74.10000 4.59000)
3	2021-01-21	6953823661	4	53.502930	7.201172	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.09500 4.58500)
4	2021-01-21	6952553476	4	26.608398	6.316406	0	ANTONIO NARIÑO	January	Thursday	POINT (-74.12100 4.59300)

OP3 = points_late[points_late['N_OPERADOR']==3][['FECHA','CLASS_LATE', 'MES','geometry']]

	FECHA	CLASS_LATE	MES	geometry
2068	2021-01-21	0	January	POINT (-74.10200 4.59200)
2073	2021-01-19	0	January	POINT (-74.10100 4.58900)
2077	2021-01-18	0	January	POINT (-74.09800 4.58000)
2079	2021-01-20	0	January	POINT (-74.10200 4.59200)
2081	2021-01-28	0	January	POINT (-74.10200 4.59200)
...	...	...	...	...
4843908	2021-01-30	0	January	POINT (-74.12600 4.47200)
4843909	2021-01-31	0	January	POINT (-74.12600 4.47200)
4843910	2021-01-31	0	January	POINT (-74.12600 4.47200)
4843911	2021-01-31	0	January	POINT (-74.12600 4.47200)
4843912	2021-01-31	0	January	POINT (-74.12600 4.47200)

41095 rows × 4 columns

Crearemos un hotplot mensual para el operador 3 y de esta forma tener un analisis por operador.

token = 'pk.eyJ1IjoiZnJhbmRhazIiLCJhIjoiY2t0cHpubW1jMDlmdzJvbXhjZ2o0eHc0ZiJ9.WiHLSWVwG9VYcmBLls4oDg' ## token para usar plotly
px.set_mapbox_access_token(token)#acceder con token

## acomodamos los datos segun las necesidades de la funciones de folium, ya que recibe una lista de cordenadas por fecha
fecha = Comparendos[['LATITUD','LONGITUD','CANTIDAD','FECHA_HORA']]
fecha = fecha.groupby('FECHA_HORA')
array=[]
time_index=[]
for key, group in fecha:
    time_index.append(str(key))
    array.append(np.array(group[['LATITUD','LONGITUD','CANTIDAD']]))
array = [i.tolist() for i in array]

## creamos un mapa con folium que sea dinamico para ver el hotsplot por dia de nuestros datos.
m = folium.Map([4.6, -74.0], tiles="stamentoner", zoom_start=12)# mapa base

hm = plugins.HeatMapWithTime(array, index=time_index, auto_play=True, max_opacity=0.3) #heatmap diario y con auto play

hm.add_to(m)# lo agregamos al mapa base

m # lo publicamos o mostramos

Make this Notebook Trusted to load map: File -> Trust Notebook

estos mapas son dinamicos y nos ayudan a mostrar de forma rapida que puede suceder con nuestros datos, pero que ocurre cuand deseamos operar con esta informacion? pues para eso tenemos muchas herramientas como usar OGR, GDAL o QGIS puro para crear esto mapas de calor o densidad

import geoplot as gplt
import geoplot.crs as gcrs

localidades = gpd.read_file('data\localidades\localidades_data.shp') #leemos el shape de localiades enriquecido

# creamos un heatmap dede geoplot
OP3_J = OP3[OP3['MES']=='January']

ax = gplt.polyplot(localidades,zorder=1)
gplt.kdeplot(OP3_J, cmap='Reds',shade=True, clip=localidades, ax=ax)

<AxesSubplot:>

# creamos un heatmap dede geoplot
OP3_S = OP3[OP3['MES']=='September']

ax = gplt.polyplot(localidades,zorder=1)
gplt.kdeplot(OP3_S, cmap='Reds',shade=True, clip=localidades, ax=ax)

<AxesSubplot:>

con estos mapas podriamos hacer diversas operaciones como sumarlos para ver frecuencias en el tiempo todo esto lo podemos hacer a una escla diaria, semana, mensual.

ahora hallaremos los cluster a nivel mensual, utilizando Kmeans, nos ayudara a clasificar cada punto dentro un grupo para ver si hay patrones espaciales

MES = OP3[['MES', 'geometry']]
MES = MES.groupby('MES')
for key, group in MES:
    gdf = MES.get_group(key)
    gdf['MES'] = gdf['MES'].astype(str) 
    a=pd.Series(gdf['geometry'].apply(lambda p: p.x))
    b=pd.Series(gdf['geometry'].apply(lambda p: p.y))
    X=np.column_stack((a,b))

X

array([[-74.098,   4.592],
       [-74.094,   4.589],
       [-74.094,   4.593],
       ...,
       [-74.126,   4.472],
       [-74.126,   4.472],
       [-74.126,   4.472]])

Usaré el método del ELBOW para determinar el número óptimo de grupos.

from sklearn.cluster import MiniBatchKMeans as KMeans
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()

def train_elbow_cluster(train,start=1,finish=10,each=1):
    wcss = {} #dicionario vacio para guardar los errores
    train = scaler.fit_transform(train)
    for i in range(start,finish,each):
        kmeans = KMeans(n_clusters = i, init = 'k-means++', random_state = 42)
        kmeans.fit(train)
        wcss[i] = kmeans.inertia_ #suma de distancias cuadradas a los centros del grupos más cercanos
    plt.title('Metodo de Elbow ')
    plt.xlabel('k= n grupos')
    plt.ylabel('WCSS')
    sns.pointplot(x=list(wcss.keys()), y=list(wcss.values()))
    plt.show()
    return wcss

train_elbow_cluster(X,finish=150,each=10)

{1: 90227.04206798342,
 11: 6066.826485932875,
 21: 2452.8866630684834,
 31: 1347.6473883638705,
 41: 1111.935567678274,
 51: 677.9085363349694,
 61: 656.1011746010108,
 71: 500.20839002660466,
 81: 522.8608182843452,
 91: 397.6404006346729,
 101: 350.5697539419241,
 111: 436.1014250045261,
 121: 303.5327723824021,
 131: 291.2312834750681,
 141: 285.3736253617679}

train_elbow_cluster(X,finish=20,each=1)

{1: 90227.04206798342,
 2: 52479.08456226006,
 3: 33051.75894019381,
 4: 17419.321418686104,
 5: 14821.24250145048,
 6: 11449.451611247408,
 7: 10783.110284358494,
 8: 8675.147416981168,
 9: 8432.834281278883,
 10: 7343.653444691406,
 11: 6066.826485932875,
 12: 5426.048165769307,
 13: 5025.335847984384,
 14: 4303.0931026706785,
 15: 4162.092746454154,
 16: 3482.384474462585,
 17: 3857.5928249627023,
 18: 3039.3769389585327,
 19: 2690.80438388283}

kmeans = KMeans(n_clusters = 16, init = 'k-means++', random_state = 42,  max_iter=400)
X = scaler.fit_transform(X)
y_kmeans = kmeans.fit_predict(X)
gdf['cluster']=y_kmeans

gdf.to_crs(epsg=4326, inplace=True)

ax = gplt.polyplot(
    localidades,
    zorder=-1,
    linewidth=1,
    projection=gcrs.Mollweide(),
    edgecolor='white',
    facecolor='lightgray',
    figsize=(12, 12)
)
gplt.pointplot(
    gdf,
    hue='cluster',
    legend=True, legend_var='hue',
    ax=ax
)

plt.title("cluster del mes")

Text(0.5, 1.0, 'cluster del mes 2010')

este mapa nos muestra que el operador 3 en septiembre tiene al rededor de 16 clusters. Todos estos analisis son replicaples para cada operdor, mes, dia y para diferentes escalas espaciales.

gdf

	MES	geometry	cluster
11848	September	POINT (-74.09800 4.59200)	15
11929	September	POINT (-74.09400 4.58900)	0
11930	September	POINT (-74.09400 4.59300)	0
11931	September	POINT (-74.09000 4.59200)	0
11933	September	POINT (-74.09400 4.58900)	0
...	...	...	...
4845546	September	POINT (-74.12600 4.47200)	13
4845547	September	POINT (-74.12600 4.47200)	13
4845548	September	POINT (-74.12600 4.47200)	13
4845549	September	POINT (-74.12600 4.47200)	13
4845550	September	POINT (-74.12600 4.47200)	13

45108 rows × 3 columns